home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / RAN3.DEM < prev    next >
Text File  |  1991-05-01  |  2KB  |  58 lines
  1. PROGRAM d7r4 (input,output);
  2. (* driver for routine RAN3 *)
  3. (* calculates pi statistically using volume of unit n-sphere *)
  4. CONST
  5.    pi=3.1415926;
  6. VAR
  7.    i,idum,j,k,jpower : integer;
  8.    x1,x2,x3,x4 : real;
  9.    iy : ARRAY [1..3] OF integer;
  10.    yprob : ARRAY [1..3] OF real;
  11.    glinext,glinextp : integer;
  12.    glma : ARRAY [1..55] OF real;
  13.  
  14. FUNCTION fnc(x1,x2,x3,x4 : real) : real;
  15. BEGIN
  16.    fnc := sqrt(sqr(x1)+sqr(x2)+sqr(x3)+sqr(x4))
  17. END;
  18.  
  19. FUNCTION twotoj(j : integer) : integer;
  20. BEGIN
  21.    IF (j=0) THEN twotoj := 1
  22.    ELSE twotoj := 2*twotoj(j-1)
  23. END;
  24.  
  25. (*$I MODFILE.PAS *)
  26. (*$I RAN3.PAS *)
  27.  
  28. BEGIN
  29.    idum := -1;
  30.    FOR i := 1 to 3 DO BEGIN
  31.       iy[i] := 0
  32.    END;
  33.    writeln;
  34.    writeln ('volume of unit n-sphere, n := 2,3,4');
  35.    writeln ('# points','     pi     ','  (4/3)*pi  ',' (1/2)*pi^2 ');
  36.    writeln;
  37.    FOR j := 1 to 13 DO BEGIN
  38.       FOR k := twotoj(j-1) to twotoj(j) DO BEGIN
  39.          x1 := ran3(idum);
  40.          x2 := ran3(idum);
  41.          x3 := ran3(idum);
  42.          x4 := ran3(idum);
  43.          IF (fnc(x1,x2,0.0,0.0) < 1.0) THEN  iy[1] := iy[1]+1;
  44.          IF (fnc(x1,x2,x3,0.0) < 1.0) THEN  iy[2] := iy[2]+1;
  45.          IF (fnc(x1,x2,x3,x4) < 1.0) THEN  iy[3] := iy[3]+1
  46.       END;
  47.       FOR i := 1 to 3 DO BEGIN
  48.          jpower := twotoj(j);
  49.          yprob[1] := 4.0*iy[1]/jpower;
  50.          yprob[2] := 8.0*iy[2]/jpower;
  51.          yprob[3] := 16.0*iy[3]/jpower
  52.       END;
  53.       writeln (jpower:6,yprob[1]:12:6,yprob[2]:12:6,yprob[3]:12:6)
  54.    END;
  55.    writeln;
  56.    writeln ('actual',pi:12:6,(4.0*pi/3.0):12:6,(0.5*sqr(pi)):12:6)
  57. END.
  58.